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We study the generation of vorticity and velocity dispersion by orbit crossing using cosmological 
numerical simulations, and calculate the backreaction of these effects on the evolution of large-scale 
density and velocity divergence power spectra. We use Delaunay tessellations to define the velocity 
field, showing that the power spectra of velocity divergence and vorticity measured in this way 
are unbiased and have better noise properties than for standard interpolation methods that deal 
with mass weighted velocities. We show that high resolution simulations are required to recover 
the correct large-scale vorticity power spectrum, while poor resolution can spuriously amplify its 
amplitude by more than one order of magnitude. We measure the scalar and vector modes of the 
stress tensor induced by orbit crossing using an adaptive technique, showing that its vector modes 
lead, when input into the vorticity evolution equation, to the same vorticity power spectrum obtained 
from the Delaunay method. We incorporate orbit crossing corrections to the evolution of large scale 
density and velocity fields in perturbation theory by using the measured stress tensor modes. We 
find that at large scales (k ~ O.f /iMpc" 1 ) vector modes have very little effect in the density power 
spectrum, while scalar modes (velocity dispersion) can induce percent level corrections at z = 0, 
particularly in the velocity divergence power spectrum. In addition, we show that the velocity power 
spectrum is smaller than predicted by linear theory until well into the nonlinear regime, with little 
contribution from virial velocities. 



I. INTRODUCTION 

The evolution of cosmological perturbations is de- 
termined, at scales larger than those where bary- 
onic physics becomes important, by the gravitational 
clustering of cold dark matter, which can be taken as 
collisionless to a very good approximation. There- 
fore, in this regime the Vlasov equation, i.e. the col- 
lisionless limit of the Boltzmann equation, describes 
the dynamics of cosmological perturbations [l| ■ 

At large scales, where orbit crossing may be ne- 
glected, the Vlasov equation reduces to the dynamics 
of a pressureless perfect fluid (hereafter PPF). The 
PPF approximation has been used extensively in an- 
alytic approaches such as standard perturbation the- 
ory (hereafter PT; see [|J for a review) and the more 
recent renormalized perturbation theory [3| (here- 
after RPT) and related techniques @, i, i, W, & 9 • 

At small scales, as the first nonlinear structures 
are formed, orbit crossing generates a nontrivial 
stress tensor (the second cumulant of the phase space 
distribution function), which leads to velocity dis- 
persion and vorticity in the dark matter distribu- 
tion. Both of these effects would not be significantly 
present otherwise; vorticity corresponds to vector 



modes which are not produced primordially (at least 
in the simplest models of inflation) and even if they 
were they decay due to the expansion of the universe, 
velocity dispersion does get generated primordially 
e.g. during thermal equilibrium of dark matter in 
the early universe but for cold dark matter candi- 
dates typical values are vanishingly small (~ 10~ 6 
km/s for WIMPs and at most ~ 10 -10 km/s for 
axions [to)) compared to typical velocity flows gen- 
erated during structure formation. 

Although a number of works have attempted to go 
beyond the PPF approximation [lj El 111, EI El 
[la E3 > there has been no quantitative estimate in 
the literature at what scale corrections to the PPF 
become important. This is the main goal of this pa- 
per. In order to achieve this, one has to compare 
PPF solutions with solutions to the Vlasov equa- 
tion. N-body simulations of collisionless cold dark 
matter attempt to solve the latter by discretizing 
the distribution function using particles that follow 
the characteristics of the Vlasov equation [65] . The 
N-body solution will differ from the PPF in regions 
where particle orbits cross, also known as "caustics" 
or "shell-crossing" in the context of the spherical col- 
lapse. This generates a nontrivial stress tensor and 



higher-order cumulants of the distribution function 
in the dark matter. We use the N-body solution to 
measure the induced stress tensor generated by orbit 
crossing and calculate from it the corrections to the 
PPF predictions for the density and velocity diver- 
gence power spectra at large scales. Recent work on 
orbit crossing has concentrated on enhancement of 
dark matter annihilation in caustics [H, [l^, [2(| ■ We 
are instead interested on the impact of orbit cross- 
ing on the large-scale dynamics, outside dark matter 
halos. 

Along the way we provide a number of results re- 
garding the nonlinear evolution of peculiar veloci- 
ties, which compared to the density field has not 
been studied in as much depth. The generation of 
vorticity and velocity dispersion impacts the recon- 
struction of primordial fluctuations from peculiar ve- 
locities [22| which assume a cold (single-stream) 
potential flow, as do other reconstruction methods 
based on galaxy positions [H, [24], [H[ . The under- 
standing of the nonlinear evolution of the volume 
weighted (as opposed to mass or galaxy weighted) 
velocity field is important, since the velocity differ- 
ence PDF is one of the building blocks that con- 
tributes to the redshift space galaxy power spec- 
trum, independent of g alaxy bias and calculable 
from first principles ( [26| ] see also [27j for recent dis- 
cussion). The study of the peculiar velocity field 
has been recently highlighted as a means to con- 
strain dark energy and large-distance modifications 
of gravity [H, EH H3|- Finally, as old and new 
techniques to measure the peculiar velocity power 
spectrum improve, some of the issues we study here 
should be important for making predictions that 
model nonlinear effects accurately for future obser- 
vations 

This paper is organized as follows. In section ITT1 
we present results for the power spectrum of ve- 
locity divergence and vorticity that follow from ap- 
plying the Delaunay method to our N-body simu- 
lations. We discuss how vorticity and velocity dis- 
persion get generated by orbit crossing in section ITTTI 
where we also propose an estimator of the stress ten- 
sor induced by orbit crossing based on an adaptive 
method. In section ITVl we extend PT to include vor- 
ticity and velocity dispersion. Finally we present our 
conclusions in section [V] 
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FIG. 1: Velocity divergence power spectrum at z — 
from 50 realizations of the LR1280 simulations. The 
symbols with error bars denote the velocity divergence 
power spectrum measured with the Delaunay method, 
normalized as in Eq. (pTJ . The blue solid line is the RPT 
prediction, and the dotted line is the linear power spec- 
trum. 



II. DIVERGENCE AND VORTICITY 
POWER SPECTRA 

A. Spatial Distribution in N-Body Simulations 

In Appendix [5] we discuss how to estimate the 
velocity field from Delaunay tessellations, also com- 
paring to more standard interpolation methods that 
deal with mass weighted velocities. We refer the 
reader to the Appendix for technical details. Here we 
show the results of applying the Delaunay method 
to estimate velocity statistics from the cosmologi- 
cal simulations described in Table fl] Note that for 
all plots in this paper, we normalize the divergence 
(9 = V • u) and the vorticity (w = V x u) so that 
they refer to the dimensionless quantities, i.e. 

V • u/Hf and V x u/Hf, (1) 

where TL — dhia/dr and / = d In D + jd In a is the 
logarithmic derivative of the linear growth factor D + 
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FIG. 2: Illustration of overdensity, divergence and vorticity in a l/i -1 Mpc thick cross-section of the simulation box 
at z = 0. The divergence and vorticity components panels correspond to the dimensionless quantities V ■ u/Tif and 
V x u/H.f. The panel labeled "Halo particle overdensity" shows the overdensity of particles belonging to dark matter 
halos with mass m > 2.4 x 10 10 /i _1 M e . 



with respect to the scale factor a. This change of 
units is convenient since in linear theory, the di- 
vergence normalized in this way equals minus the 
dimensionless overdensity, i.e. 8 = —TLfS, with 
S = Sp/p. 

Figure [1] shows the average, over the 50 realiza- 
tions of run LR1280, of the divergence power spec- 
trum at z — compared with the 2-loop RPT pre- 
diction. The divergence power spectrum behaves as 
expected theoretically, with suppressed growth com- 
pared to linear theory, although one can notice sig- 
nificant deviations from RPT, which will be explored 
in detail elsewhere. Note that the non-linear effects 
in the power spectrum are observable on scales with 
k > 0.01 hr 1 Mpc, unlike the case for the density 
field. This is expected for two reasons: first, the 
velocity divergence propagator decays faster than 
for the density field, damping the linear spectrum 
faster with k [481 ] : second, the mode coupling power 
generated at small scales is smaller than for the 
density field, avoiding the accidental cancellation of 



nonlinear effects present in the density power spec- 
trum [49( . This qualitative behavior is also predicted 
by standard PT [26| . Clearly, as discussed in [2(| , 
assuming that density and velocity divergence power 
spectra are equal (as often done for redshift distor- 
tions) is not a very good approximation, even at 
large scales. 

Figure [5] illustrates the spatial distribution of ve- 
locity and density estimations from the HR160 run 
at z = 0. The panels show overdensity, divergence 
and vorticity on a 1 h^ 1 Mpc-thick cross-section of 
the simulation. Also, the overdensity corresponding 
to halo particles (particles inside dark matter halos) 
is shown. Even though the overdensity field can take 
on values up to a few hundreds, its scale was chosen 
to go up to 6 = 3 because the dark matter halos 
are small compared to the scale of this figure and 
increasing the upper scale limit would just hide the 
lower density structures. 

We can see that the divergence field is, not surpris- 
ingly, remarkably similar to the density field. How- 
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B. Dependence on Mass Resolution 
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FIG. 3: Dependence on mass resolution of the velocity 
divergence and vorticity power spectra. Particle masses 
are labeled in units of 10 10 h -1 Mq, see TableUfor more 
details on the simulations. While the divergence power 
spectrum does not depend on mass resolution, the vor- 
ticity power spectrum does show significant sensitivity. 
However, as mass resolution increases it converges (when 
m par is below ~ 10 9 /i. -1 Mq) to a stable answer. 



ever, the structures in the velocity divergence have 
in general lower amplitude and are more extended in 
space, as expected from the power spectrum results 
discussed above. It is interesting to note that, at the 
halo positions, the divergence tends to be smaller 
than in the still collapsing regions, as it should be. 
On the other hand, the vorticity field fluctuates in 
sign on scales of the order of ~ 1 h^ 1 Mpc (roughly 
as expected from theoretical arguments, see Fig. 8 
in (501), and it is concentrated on collapsing regions, 
where shell-crossing is currently occurring. There 
are no large-scale coherent fluctuations in vorticity, 
so we expect the vorticity power spectrum to be 
much smaller than the divergence power on large 
scales, as we now discuss. 



Figure [3] shows the power spectrum of divergence 
and vorticity obtained from the Delaunay method 
from different simulations (see Table Hj). The veloc- 
ity field is dominated, especially on large scales, by 
its irrotational component, consistent with the spa- 
tial distribution seen in Fig. [2] We see from Fig. [3] 
that the divergence power spectra measured over a 
broad range of volume, number of particles and mass 
resolution simulations match consistently. 

The estimate of the vorticity power spectrum, on 
the other hand, appears not to be so robust: it 
shows a clear monotonic dependence on the mass 
resolution. We verified that this dependence was 
not an artifact of the Delaunay method by com- 
paring these results to the ones obtained from the 
CIC mass-weighted scheme. We observed that these 
two methods agree on the mass resolution depen- 
dence of the vorticity power spectrum (not shown in 
Fig. [3] for clarity). Thus we believe the dependence 
on mass resolution of the measured vorticity is real 
and may be due to insufficient sampling of collaps- 
ing regions 66]. However, as the particle mass goes 
below 10 9 h 1 M , the vorticity power spec- 

trum eventually converges. 

Also, we check for aliasing effects, discussed in de- 
tail in the Appendix. Our estimates for the spurious 
aliased vorticity based on Eq. (|A13|) are at least two 
orders of magnitude lower than the measured vortic- 
ity from the simulations. Although this shows the 
spurious vorticity it is not a sampling issue in the 
measurement of the power spectrum, aliasing may 
be an issue in the low resolution simulations (which 
have a coarser PM grid) during time evolution, since 
the spurious power is close to the expected k 2 behav- 
ior (see section I A 41 in the Appendix) at all scales. 
From here on, our results will be based only on the 
higher mass resolution runs. 

Regarding the dependence on resolution of the di- 
vergence power spectrum, close inspection of Fig. |3] 
seems to indicate that the higher mass resolution 
runs have a larger power spectrum than lower mass 
resolution runs by about 5 — 10%. However, one 
must keep in mind the higher mass resolution runs 
have substantially smaller box sizes (see Table HJ. 
For box sizes smaller than about ~ 300 hr 1 Mpc, 
RPT predicts that the propagator is seriously af- 
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Name 




■^V; part 




N realizations 


softening 


LR1280 


1280 


640 3 


59.94 


50 


0.07 


LR512 


512 


256 3 


59.94 


1 


0.2 


MR512 


512 


512 3 


7.49 


1 


0.04 


MR320 


320 


640 3 


0.94 


1 


0.015 


HR160 


160 


640 3 


0.12 


1 


0.00625 


SHR160 


160 


1024 3 


0.029 


1 


0.00625 



TABLE I: All our simulations have fi m = 0.27, Q A = 0.73, Q b = 0.046, h = 0.72 and a 8 (z = 0) = 0.9. They were run 
using the Gadget2 code [43]. L bo x is in units of h' 1 Mpc and m par is the particle mass in units of 10 1C 'ft. -1 Mq. 



fected by the finite volume of the simulation (see 
Fig. 6 in [HI). For such boxes, the damping of the 
linear spectrum by the propagator is much less se- 
vere, and while the mode coupling power is some- 
what larger, it cannot compete with the nearly ex- 
ponential scale dependence of the propagator. Thus 
one expects to see higher divergence power in smaller 
boxes. This is confirmed further by looking at simu- 
lations of same box size (LR512 and MR512, on one 
hand, and HR160 and SHR160 on the other). The 
ratio of the power spectra in boxes of the same size 
but very different mass resolutions does not show 
any significant (percent-level) deviation from unity. 



Finally, note that finite volume effects are not ex- 
pected to affect the vorticity power spectrum, as it is 
dominated by small scale structures. That the vor- 
ticity is sensitive to mass resolution rather than box 
size is clear from comparing the LR512 and MR512 
results, which differ by a factor of 8 in mass reso- 
lution but have the same box size. Figure [3J shows 
their vorticity power spectra differ by a factor of 
about four. 



It is worth noting that the velocity power spec- 
trum obeys P v (k) = k 2 [P e (k) + P w (k)}. Thus, the 
resolution dependence of P w seen in Fig. [3] means 
that when P w is comparable to Pg, a similar depen- 
dence on resolution affects the velocity power spec- 
trum. Therefore, spurious vorticity can lead to an 
overestimate of the velocity power spectrum at non- 
linear scales. 
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FIG. 4: Time dependence of the divergence and vortic- 
ity power spectra. The divergence power spectrum at 
z — 1 and z — 3 are linearly extrapolated to z = for 
comparison. The vorticity power spectrum was similarly 
scaled using Eq. (|3]) with n w = 7. In the non-linear 
regime, both divergence and vorticity grow slower than 
the large-scale extrapolation. 



C. Time Dependence 

As will be shown later, in order to calculate how 
much vorticity affects the evolution of the density 
power spectrum, it is necessary to determine the 
time dependence of the vorticity power spectrum. In 
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linear theory, the divergence power spectrum evolves 
according to 



P e (k,z) ~ [D + {z)f P eo (k), 



(2) 



where Peo(k) is the initial (post-recombination) di- 
vergence power spectrum and D+[z) is the linear 
growth factor measured away from the initial condi- 
tions. For the vorticity power spectrum, we propose 
that in the large-scale limit 



P w {k,z) oc [D + (z)Y 



(3) 



From the vorticity power spectra estimates at 
z = 0, 1,3, we find that the best fit for Eq. is 
given by n w — 7±0.3. Figured] shows the divergence 
power spectrum at z = 0,1,3 extrapolated to z = 
using Eq. @, and the vorticity power spectrum ex- 
trapolated using Eq. Q with n w — 7. It can be seen 
that at large scales the extrapolated outputs agree 
very well while in the non-linear regime these simple 
scalings break down, as expected, with the growth 
slowing down compared to large scales. In the case 
of the divergence power spectrum, this behavior can 
be understood reasonably well from RPT, see Fig.[TJ 
A detailed discussion of velocity statistics and RPT 
will be presented elsewhere. 

For the vorticity, little is known from first prin- 
ciples. The exception is the work in [50| , were the 
rms vorticity is calculated for CDM spectra during 
first orbit crossing using the Zel'dovich approxima- 
tion (see also [13[ for a calculation of rms vorticty 
for scale- free models). They found their estimates 
were much smaller than found in simulations, which 
given their mass resolution at the time is not sur- 
prising (see Fig. (31). However, one can very roughly 
estimate the vorticity power spectrum at large scales 
generated by orbit crossing by postulating that orbit 
mixing creates in such regions a velocity field that is 
approximately the result of mass weighting the single 
stream velocities (see discussion related to Eq.B5lbe- 
low). Then we expect w ~ (1 + <5)~ 1 / t ,V x [(l + d)v], 
where /„ is only nonzero in regions where orbit cross- 
ing occurs, and on average can be thought as the 
fraction of the volume that undergoes orbit crossing, 
an increasing function of time. Then the vorticity 
power spectrum reads [46| 



P w (k,z) ~ [f v {z)f I d 3 q 



(k x q) 2 



P 5 (\k-q\)Pg(q) 



(k - q) 2 



P x (|k-q|)P x (g) 



(4) 



where P x is the cross spectrum between 5 and 9. In 
the low-fc limit this reduces to 



P w {k,z) ~ [f v (z)} 2 / d 3 q 



(k x q) 2 



P&(q)Pe(q) 



P x (q)P x (q) cx k 2 [f v (z)} 2 [D + (z)} 6 , 



(5) 



where in the last step we have assumed the veloci- 
ties are normalized as in Eq. ((T|) and used that the 
square brackets vanish in linear theory, so the lead- 
ing nonzero contribution comes from one-loop PT 
(which induces a time dependence in the power 
spectra beyond leading order). Despite the crude 
approximations made in arriving to Eq. ((5]), the 
scale and time dependence of the large-scale vortic- 
ity power spectrum seen in Fig. [4] may be explained 
qualitatively along these lines. 



D. Impact of Virial Velocities 

In Figure [21 it can be observed that the velocity 
field is rotational in high-density collapsing regions. 
Compare, for instance, the lower panels against the 
top right panel for which only particles belonging to 
halos are shown. On the upper left corner of the 
simulation box, there is a large filamentary struc- 
ture. We can see that the vorticity occurs mainly on 
the outskirts of virialized objects. This suggests that 
the fraction of the vorticity power spectrum com- 
ing from virialized regions themselves is not very 
big. To check this, we took the HR160 simula- 
tion and replaced the particle velocities belonging 
to halos by the center-of-mass velocity of the cor- 
responding halo, thus eliminating the velocity dis- 
persion of all halos. We measured divergence and 
vorticity power spectra and compared them to those 
of the unmodified HR160 simulation. The results 
are shown in Figure [5l It can be seen that at the 
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FIG. 5: Comparison between divergence and vorticity 
power spectra of the simulation to the power spectra ob- 
tained by replacing velocities of particles inside halos by 
the center of mass velocity of the parent halo, thus set- 
ting virial velocities to zero. We can see that the diver- 
gence is mostly unaffected, while the vorticity differences 
are less than 5%. 



scales we probe the divergence power spectrum is 
essentially not affected by the virial velocities, and 
the vorticity power spectrum is reduced by less than 
5%. It is important though to keep in mind that our 
measurements on the HR160 simulations are done 
in a grid of size 160, so the contribution from scales 
less than ~ 1 hr 1 Mpc are not included; this roughly 
corresponds to ignoring halos of m < 10 14 /i~ 1 M Q . 

That the vorticity power spectrum is not very sen- 
sitive to virial velocities may be understood by con- 
sidering the vorticity evolution equation (0], see also 
Eq. [37J below) 

<9w fl \ 

— +Hw- V x [u x w] = V x -V • (pa) , (6) 

OT \P J 

where pa^ is the stress tensor induced by orbit cross- 
ing. In a virialized object, where the phase-space 
distribution function is approximately Maxwellian 
with velocity dispersion related to halo mass through 
a 2 ir ~ Gm/r V i r (r V i r oc m 1 / 3 ), we expect the stress 



tensor to be reasonably well approximated by an 
equation of state of the form poij ~ —pdij, where 
p is a density-dependent pressure; in fact, for an 
isothermal sphere p oc pa 2 ir with a 2 ir independent 
of spatial position. In that case, the forcing term of 
Eq. ((6|) can be written as 



V x -V • (pa) 



x Vp; 

p 2 



0. 



(7) 



where the last step follows from the density- 
dependence of the pressure. Therefore, although 
these approximations are not totally realistic in 
practice, they may help explain small vorticity 
sourcing from virialized dark matter halos at the 
scales we probe. 



III. GENERATION OF VORTICITY AND 
VELOCITY DISPERSION 

A. Beyond PPF 

One of the main goals of this paper is to esti- 
mate the corrections to the pressureless perfect fluid 
(PPF) predictions for power spectra of density and 
velocity fields at large scales due to orbit crossing at 
small scales. In order to do so, we have to estimate 
the corrections to the equations of motion beyond 
PPF that result from solving the Vlasov equation 
for the phase-space distribution function (hereafter 
DF) /(x,p,r), 



df p df 
+ _ . y/ - aV(j> ■ — 
or a op 



(8) 



where p is the momentum per unit mass and the 
gravitational potential. Equation ([8]) says that the 
DF is conserved (df/dr = 0) along its characteris- 
tics, 



fix 
dr 



P 

a 



— = -aS7(j>, 



(9) 



which are the Hamilton equations of motion, that 
can be combined to give the familiar result, 



dr 2 dr 



-V<t>. 



(10) 
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Cosmological N-body simulations solve the Vlasov 
equation by discretizing the DF using particles that 
follow the characteristics, Eqs. ^ or (fTU|) . To make 
connection with the PPF equations of motion, one 
may take moments (or, more precisely, cumulants) 
of the Vlasov equation. The first few cumulants of 
the DF are the (comoving) density field, 



(1 + 8) 



Hp) 



d 3 P , 



(11) 



where to simplify notation we avoid displaying the 
space and time arguments everywhere; the velocity 
field u, 

(l + 5)u = J f(p)^d 3 p, (12) 

and the stress tensor Ty = (1 + 5)<Jij, 

(1 + 5) oij = J /(p) ^ d 3 P -(1 + 5) u iUj , (13) 

where the velocity dispersion tensor cry describes 
isotropic and anisotropic velocity dispersion. Be- 
fore we derive equations of motions for these quanti- 
ties, it is useful to introduce the cumulant generat- 
ing function, which generates all these objects. As 
it is usual in statistics of large-scale structure (see 
e.g. 0), the cumulant generating function C is given 
in terms of the moment generating function M. by 

C(l) = lnX(l), M(\) = J e^/ /(p) d 3 p, (14) 

where moments are obtained by successive deriva- 
tives of M. with respect to the external parameter 
I, 

(V kl ...V hn M) = (l + 6) m^ An (15) 



where (. . .)o means evaluating quantities at 1 = 0, 

and mW = 1, M = (1 + 8), mf ] = u u m\f = 
U{Uj + Uij. Cumulants are statistically independent 
objects at each order and can be obtained similarly 
by differentiation of C, 



(V iii ...V; in C) = c^ n , 



(16) 



„(2) 



with c (0) = C = ln(l + 5), c\ ' = m, c^' = a i} . 



From the Vlasov equation, Eq. (jHJ), and Eq. (|T3]1 
it is straightforward to derive equations of motion 
for the generating functions. For C we have, 



dC 



n(i- voc + vc • ViC + (v • voc 



-1- V0. 

(17) 

This is a nonlinear partial differential equation for 
C(x, r, 1); however, all we are interested in is what 
happens in the neighborhood of 1 = 0, i.e. the 
derivatives of C at 1 = 0, see Eq. (|16[) . 

The equations of motion beyond the PPF approx- 
imation readily follow from Eq. (fT"7)) . Setting 1 = 
we obtain the continuity equation, 



g + V-[(l + 5)u] = 0, 



(18) 



whereas taking the first derivative we obtain mo- 
mentum conservation, 



duj 
dr 



TLui + (u • V)iti 



-W<j>--V j (pa ij ), (19) 



where p = 1 + 5, while the evolution of the veloc- 
ity dispersion tensor is obtained from Eq. (|17p by 
applying second derivatives, 



d<7jj 



2'Hcr ij + (u • V)cr. 



(20) 



1. 



CTjk^kUi + CTik^kUj = --Vfc(pITj fc ), 

- J 3 ) 



c\s k is the third cumulant of the DF, 



where ri^fc 

see Eq. (fTo]) . By applying successive derivatives with 
respect to 1 in Eq. (jTTJ) one thus generates an infinite 
hierarchy of equations of motion for the cumulants of 
the DF (hereafter cumulant hierarchy). The hierar- 
chy is infinite because at finite order is never closed, 
the cumulant of order n depends on that of order 
n + 1. 



B. The Cumulant Hierarchy and Orbit 
Crossing 

Such an infinite hierarchy is very difficult to solve. 
The PPF approximation truncates the hierarchy as- 
suming that the second and higher order cumulants 
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of the DF are zero, thus <7y = 0, Tlijk — in 
Eqs. (|19H20p and so on. This is equivalent to as- 
suming the DF takes the form, 

/(x, p, r) = [1 + 5(x, t)] S D [p-a u(x, r)], (21) 

for which C(l) = ln(l+#)+l-u, and clearly all cumu- 
lants of order larger than one vanish. Note that the 
PPF approximation appears to be self- consistent, 
i.e. assuming that aij and higher order cumulants 
vanish at a given time is preserved by the hierarchy. 
This is so because there are no linear or nonlinear 
terms in the equations of motion for such cumulants 
that solely involve the density and/or velocity fields 
as sources. This can be readily seen from the struc- 
ture of Eq. (117|) . after operating with two or more 
derivatives Vi. In other words, if C initially only 
contains linear terms in 1, Eq. (|17p will not generate 
higher powers in 1. 

This situation is, however, unstable under pertur- 
bations. If somehow C develops a quadratic con- 
tribution in 1, then the nonlinear term in Eq. (|17p 
generates a cubic term, and this in turn generates 
higher orders, and so on. Therefore, once velocity 
dispersion "turns on" , all higher order cumulants do 
so as well. Thus a priori it is not a self-consistent 
truncation to include a non-zero o~ij and ignore Fl^-fc 
(which is sourced by terms solely dependent on ) 
and higher-order cumulants. This truncation may, 
however, become a good approximation in some sit- 
uations, e.g. at large scales. 

Physically it is expected that even for perfectly 
"cold" initial conditions where the DF is given ini- 
tially by Eq. (|21[) . orbit crossing during time evo- 
lution will generate a nontrivial stress tensor and 
higher-order cumulants, while as discussed above the 
cumulant hierarchy does not seem to allow for this. 
Given that such a result from the cumulant hier- 
archy is unstable to small perturbations away from 
cold initial conditions, any subtlety in going from 
the Vlasov equation to Eq. (flT)) may alter the con- 
clusions. A more careful look at orbit crossing in this 
context shows that this suspicion is well founded. 

To see how orbit crossing generates a nontrivial 
DF from cold initial conditions, consider the formal 
solution of the Vlasov equation expressing the con- 
servation of the DF along the characteristics, 

/(x,p,r) =/„(Xo,Po), (22) 
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FIG. 6: Phase space sketch of generation of multiple 
streams due to orbit crossing. The top panel shows the 
zero width DF after evolution from cold initial conditions 
before orbit crossing. The bottom panel shows the DF 
after orbit crossing (which occurs at the moment when 
"stream 2" is infinitesimal and perpendicular to the x- 
axis, while the two vertical lines coincide, as well as the 
filled circles). In between the two vertical lines there 
are three "streams": that's the region of space where 
multistreaming is present, and e.g. velocity dispersion is 
generated. The intersections of the vertical lines with the 
streams at the filled circles denote where the derivative 
of the curve is infinite and thus at such positions the 
density field (projection of the DF onto the x-axis) is 
singular. 

where /o is the initial DF, and, 

X = X (x,p,r), P =P (x,p,r) (23) 

are the initial positions and momenta which when 
evolved by the equations of motion until time r 
(Eqs. O lead to x and p. That is, time evolution 
maps (Xq,Po) to (x, p) at time t, and this map- 
ping is invertible because in phase space trajectories 
never intersect for a Hamiltonian flow. 

However, orbits clearly can (and do) cross in con- 
figuration space, i.e. at time r and position x there 
may be more than one orbit (with different p's) that 
trace back to different initial conditions (Xq,Pq)- 
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If we start from cold initial conditions, /q satisfies 
Eq. (f2Tj) . and after time evolution the DF reads, from 
Eq. dUD 



/(x, p, r) = [1 + J (Xo)] <J D [P - Uo(Xo)], (24) 

where we have set ag = 1 and So and Uo are the 
initial density and velocity fields obtained, typically, 
from Gaussian random field initial conditions. Now 
we are ready to see the effect of orbit crossing on 
the cumulant generating function C(l), Eq. (fT4")l . As 
long as orbits do not cross, since /q has zero width 
then at fixed (x, r) there is a unique p that con- 
tributes to the momentum integral (see top panel in 
Fig. [S]). Thus the argument of the delta function 
in Eq. (|24|) can be linearly related to p, and C pre- 
serves its linear dependence on 1 and no higher-order 
cumulants are generated. 

However, as soon as orbits cross, there are many 
p's at fixed (x, r) and thus many roots of the argu- 
ment of the delta function in Eq. (f!M]) . each of them 
corresponding to one "stream" , see bottom panel in 
Fig.[SJ As a result, the cumulant generating function 
reads instead 



C(x,t,1)=1h[ Yl ( 1 + < y 



(25) 



streams at x 



where we have written schematically 8 S and u s for 
the density and velocity fields of each stream, which 
can be obtained by projecting each piece of the DF 
separately, see Fig. [6j Clearly, if the number of 
streams is larger than one, C is a fully nonlinear 
function of 1 and all cumulants have been gener- 
ated simultaneously by orbit crossing. Note that 
the number of streams at position x and time r is a 
random field that depends on initial conditions and 
cosmological parameters, see e.g. [5l[ for a calcula- 
tion of the mean number of streams from Gaussian 
initial conditions. 

We now see there is an apparent contradiction be- 
tween the Vlasov formal solution and the time evo- 
lution of the cumulant hierarchy. The root of this 
can be traced back to the well known fact that at or- 
bit crossing the density field becomes singular (67| . 
therefore the cumulant hierarchy must be supple- 
mented with a regularization procedure in order to 
follow time evolution after orbit crossing. This must 



restore the agreement with the formal solution of the 
Vlasov equation, which does not suffer from develop- 
ment of singularities (which invalidate the hierarchy 
because projection in momenta does not commute 
with time evolution). 

Does this matter in practice? After all CDM 
has a very small but non-zero velocity dispersion, 
which will automatically regularize the singularities 
in the cumulant hierarchy that arise at orbit cross- 
ing. However, it seems unlikely that the final answer 
for large-scale density and velocity fields after orbit 
crossing should depend sensitively on the velocity 
dispersion of the CDM particles (if this were so CDM 
N-body simulations would be incorrect); although of 
course such effects are important for warm dark mat- 
ter candidates. Rather, it should be the self-gravity 
of regions that undergo orbit crossing that leads self- 
consistently to velocity dispersion and higher-order 
cumulants that regularizes the singularities. 

To carry out such self-consistent regularization, 
one can proceed in at least two ways: 1) introduce 
some non-zero initial velocity dispersion <jq and use 
the hierarchy, which does not develop singularities in 
this case, to evolve the system forward in time. To 
make predictions for systems with negligible initial 
velocity dispersion such as CDM one must take into 
account that the ag — > + limit is nontrivial and one 
should get finite corrections for infinitesimal Co- 2) 
since mass does not diverge in caustics, one can work 
with coarse grained variables (smoothed over some 
small scale). To find the equations of motion for 
smoothed quantities one must take into account that 
coarse graining does not commute with time evolu- 
tion, and the coarse graining scale must be picked so 
that physics at large scales is invariant with respect 
to the degrees of freedom that are integrated over at 
small scales. See [l2| for some steps in this direction. 

In either case, the net result of regularization 
is that higher-order cumulants of the DF will be 
sourced by density and velocity fields, which leads 
to an effective equation of state for dark matter. On 
the other hand, one would still have to implement a 
consistent closure of the hierarchy. 

In this paper we proceed in a different fashion, by 
measuring the stress tensor directly from numerical 
simulations. We then close the hierarchy by using 
the measured stress tensor in the momentum con- 
servation equation, Eq. (|19[) . We start our measure- 
ments at z = 3 and assume that the dark matter has 
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undergone sufficient shell crossing before z = 3 that 
future caustics are not singular, and thus Eqs. ([TBI - 
I19p are valid. To extrapolate backwards in time we 
use the time dependence found from z = 3 to z = 
and assume it is valid earlier down to the initial con- 
ditions. Since we are interested in the large-scale 
statistics of density and velocity fields, the devia- 
tions from the PPF approximation are very small 
before z — 3, thus this should be a reasonable ap- 
proach. Before we describe how we implement such 
a procedure in detail, we must explain how we esti- 
mate the stress tensor (Ty = p/Jij ) from the numer- 
ical simulations. 



C. Estimating The Stress Tensor 

In the Appendix [Al we discuss how the Delaunay 
method is a reliable algorithm for recovering the ve- 
locity field from numerical simulations. However, as 
we will now argue, it is not adequate for estimat- 
ing the velocity dispersion tensor er^ . The Delau- 
nay method is optimized for interpolating the veloc- 
ity field on arbitrary points in the simulation. This 
procedure recovers a continuous field that is to be 
interpreted as the mean velocity field. Then, if one 
tried to estimate the dispersion in a volume A V with 
an estimator of the form 

CTDcl = j ' d 3 X UDcl(x) ®u Dcl (x) - 



UDcl <8> UDol 



(26) 



where Urj e i(x) is the Delaunay-interpolated velocity 
field and ud i is the average on the volume Ay of 
UDci(x), the integral would be dominated by the De- 
launay linear approximation to the mean velocity, 
contaminating the true velocity dispersion coming 
from multi-streaming. 

Our strategy to estimate the velocity dispersion 
is based on the fact that particles in numerical sim- 
ulations are a sample of the phase-space distribu- 
tion function /(x, p, t) (see Eq.[S]). Consider, for in- 
stance, a small volume AV in which the distribution 
function is nearly constant (V/ « 0). Our ansatz 
is that the simulation particles in AV are sampled 
from the probability distribution given /av(p, t) by 



/av(p,t) 



1 



pAV 



d 3 xf(x,p,T), (27) 



where p is the density in the small volume (which is 
constant due to the ansatz). Therefore, the volume- 
averaged velocity dispersion tensor <xy in that small 
volume can be written as 



1 

AV 



Ay 



UjU j 



(28) 



where Ui and Uj are the mean velocity field in the 
volume AV (the velocity field is also constant on AV 
due to the ansatz). The term in square brackets is 
the one defined in Eq. (|2"7} . Since we assume the 
particles in AV are sampled from that distribution, 
we can write 



i N 

1 \ "* («) (") 

n=l 



UiUj, 



(29) 



AV 



where the sum is over the N particles in the small 
volume AV and u(") is the velocity of the n-th par- 
ticle. 

Equation (|29p will serve us as an estimator for 
the velocity dispersion on a small region of con- 
stant phase-space distribution function. To obtain 
the volume- weighted dispersion on a larger volume, 
we simply break that volume into regions where 
Eq. (|2"9"|) is valid and then volume-average the re- 
sults. For a different approach to estima ting velocity 
dispersion from simulations see l53l . [HJ 

In practice, we want to compute the volume- 
averaged velocity dispersion on a grid, similarly to 
Eq. (|A6|) for the mean velocity. We use the following 
recursive algorithm to find that average in a given 
cell: 

i) First determine whether cell is homogeneous 
enough. To do this compute density pi 
in each octant of the current cell. Then 
compute density standard deviation, s = 

I Elite -p) 2 - 

ii) If s < sthreshoid or N < N paltm m, we can use 
Eq. {29]). Then, compute u = ^ £n=i u| n) 
and a = 4 £n=i u<") ® - u <g> u. 
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iii) Else, the cell is not homogeneous. Following 
this three-step algorithm, calculate recursively 
the velocity dispersion crW of each octant, then 

Here, Sthreshoid is the minimum standard deviation 
of the octants density that defines the criterion for 
homogeneity (we use typically values of 30%-50% of 
the octants mean density). The purpose of -ZVp ar t m in 
is to avoid breaking cells with too few particles into 
octants (we use N paltm m — 5). Note that this al- 
gorithm is adaptive: it resolves the inhomogeneous 
dense regions to find the subregions where Eq. (f2"5)) 
holds. The method has two free parameters, which 
varied on a reasonable range produce shifts in the 
observed velocity dispersion tensor up to 40% (this 
will suffice for our purposes, as we shall see be- 
low). This is understandable, since these parameters 
compensate for the finite resolution (and sampling) 
of the phase-space distribution function. Another 
drawback of the method is that it suffers from Pois- 
son noise in low density regions. However, on large 
scales, since we are averaging over many cells, we 
expect these errors to be rather small. 

The stress tensor estimated in this way satisfies a 
nontrivial sanity check, as we will show below. Its 
vector modes source the growth of vorticity, which 
otherwise would not grow. We will compare the 
growth of the vorticity power spectrum from the vec- 
tor modes of the measured stress tensor and see that 
it agrees with direct measurements of the vorticity 
power spectrum using the Delaunay method. Before 
we show this, we need to explain how we incorporate 
the measured stress tensor into the standard calcu- 
lation of large scale evolution of density and velocity 
fields using perturbation theory. 

IV. LARGE-SCALE CORRECTIONS TO 
PPF 

A. Scalar- Vector Decomposition 

We are now ready to see the effects of orbit cross- 
ing in the large-scale evolution of density and ve- 
locity fields. We will use the equations of motion 
Eqs. (|18H19p . supplemented by the Poisson equation, 

vV = |w 2 n m 5, (30) 



and solve for the coupled system of S and u, and 
treat the stress tensor as a forcing term, with known 
scale dependence and time evolution obtained from 
measurements in the simulations. 

A non-zero stress tensor leads to two new effects 
in the evolution of large-scale structure. We can de- 
compose the velocity vector into scalar and vector 
modes, where the scalar mode is the velocity diver- 
gence (corresponding to the velocity parallel to the 
wave vector k) and the vector modes correspond to 
the vorticity (or the two components of u perpendic- 
ular to k). In the pressureless perfect fluid (PPF) 
approximation, the divergence grows (since it cou- 
ples to the gravitational potential), while any pri- 
mordial vorticity decays in linear theory due to the 
expansion of the universe, since it does not couple 
to the gravitational force because it is conservative. 
Nonlinear terms can amplify vorticity but cannot 
create vorticity if there is no primordial contribu- 
tion (as we assume in this paper). 

However, orbit crossing induces a nontrivial stress 
tensor and this will modify the evolution of the scalar 
and vector modes. From Eq. (Til?)) we see that the 
new source for scalar and vector modes is the quan- 
tity 

9i(x,r) = -Vj(po-ij). (31) 

More precisely we can decompose this into scalar 
and vector sources, 

qe = V • q, q w = V X q, (32) 

respectively. We will loosely call the correction from 
cry to the scalar modes (due to q$) "velocity disper- 
sion". Of course the velocity dispersion tensor cr^ 
affects the vector modes as well. However, in the 
simplest case of a diagonal ov,- which depends only 
on density, only qg survives. Note also that some- 
times we refer to the stress tensor = pa^ instead 
of the velocity dispersion tensor er^ for conciseness. 
Finally, for simplicity we will refer to the decomposi- 
tion of q into qg and q w as "decomposing the stress 
tensor into scalar and vector modes" . 

Let us first write the linearized version of Eqs. (fTSl - 
taking into account Eqs. (|31ti32|) . As usual, it is 
simplest to work with a different time variable, 

V = lnD + (r), (33) 
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FIG. 7: Time dependence of the power spectra of the 
scalar and vector forcing terms q e (top three lines) and 
(bottom), see Eq. ([32]) . As in Fig.H] each power spec- 
trum is linearly extrapolated to z — 0. In this case, the 
time evolution of both forcing terms is fitted by P(k, z) = 
[D + {z)/D+{0)] n - d P{k,0), with n vd = 6.5 ± 0.5. 

where D + is the linear growth factor, and scale the 
velocity and stress tensor so that, 

u->-W/u, an ^{Hffan, (34) 

and thus also qe — > (T~if) 2 qe and same for q m . Here 
/ = d\nD+/d\na w ft 5 / 9 for ACDM models. As- 
suming that f 2 w f2 m , the linearized equations of 
motion can be written after these transformations 
in the simple form, 

d v 5-e = 0, (35) 

d,0+-- Y = q e , (36) 
w 

d n w + — = q w . (37) 

Figure [7] shows the power spectra corresponding 
to the two forcing terms and their time depen- 
dence, measured from the HR160 simulation using 
the method described in section IHI CI Similarly to 
section Hi CI we fit for a time evolution of the form 

P q (k,z)<x[D + (z)] n »*, (38) 




k (h Mpc-') 



FIG. 8: Comparison between the Delaunay-estimated 
vorticity power spectrum (solid line) to the linear theory 
prediction from solving Eq. (|39|) (dotted line) for red- 
shifts z = 0, z = 1 and z = 3. 



where P q stands for both P qe and P qw . We found 
that the best fit is n vd = 6.5 ± 0.5, although the 
quality of the fit is not as good as in the case of 
the vorticity (see Fig. [4j. The reason for this may 
be shot-noise error coming from poorly sampled re- 
gions, where the error gets amplified by the 1/(1 + 6) 
factor in Eq. (f3l"j) . 

It is interesting to note that Eq. (|37|) provides 
us with a non-trivial consistency check between the 
vorticity power spectrum measured by the Delaunay 
method, and the adaptive method described in sec- 
tion UTTC] from which we measured the stress tensor 
and estimated the forcing term for vector modes. 
The vorticity and vorticity-forcing terms measured 
from the simulation should be consistent with the 
time evolution given by Eq. (f3"T| . Since this equation 
is decoupled from the other two equations (scalar 
and vector modes do not couple in linear theory) , it 
can be easily solved. Ignoring the decaying mode, 
the linear theory solution for the vorticity power 
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divergence 




scalar mode of the stress tensor corrects the PPF ap- 
proximation already at the linear level, whereas the 
vector modes are decoupled in linear theory and cor- 
rect the PPF at higher-order in PT. In this section 
we estimate the corrections due to the scalar mode 
qg (roughly speaking, velocity dispersion), while in 
the next section we tackle the corrections induced 
by q w at leading order in nonlinear PT. Since these 
deviations are small at large scales we can consider 
them separately. 

The scalar mode correction can be included by 
writing the modified linear theory of Eqs. (|35ti36|) in 
a compact form by using a two-component object 
ipi = 5, V ; 2 = that obeys the linear equations of 
motion, 

d v tf; a (k, 77) + n ab ip b (k, 77) = Q a (k, 77), (40) 
where fl a b is the 2x2 matrix, 



FIG. 9: Correction to the PPF approximation for the 
velocity divergence (three top lines) and density power 
spectrum (three bottom lines) due to velocity dispersion 
at redshifts z — (solid), z — 0.5 (dashed) and z = 1 
(dotted). Note that the actual correction is negative in 
all cases, we plot their absolute values. These correc- 
tions are computed in linear theory, Eqs. (|45[1 and (|48p . 
thus extrapolation well beyond A; ~ 0.1 /iMpc" is only 
illustrative. 



spectrum reads 



P w (k)=(—^—) P qw {k). 
\n vd + lj 



(39) 



Figure [8] shows the results of this consistency check. 
In it, we show the measured left and right hand sides 
of Eq. (f3"9"]l for redshifts z = 0, 1,3. The agreement 
in all cases is very good, improving, as expected, for 
higher redshifts. 



B. PT + Velocity Dispersion 

We are interested in estimating the large-scale cor- 
rections to the PPF approximation due to the orbit- 
crossing induced qg and q^,. As we can see from 
the linearized equations of motion, Eqs. (|35ti37[) . the 



-1 
_ 3 1 
2 2 



(41) 



and Q(k, 77) = (0, q e (k, 77)). The formal solution to 
these equations can be written as 

r 

ipai^v) = 9ab(v)<hO s -) + / H 'g a b(v~v')Qb(Kv'), 
Jo 

(42) 

where <fi represents the initial conditions and g a b is 
the linear propagator (48| . 



9ab(v) 



3 2 




(43) 



5 I 3 2 I 5 

Then, the density field in linear theory is given by 

«(M) = «M) + Kd /2^)K!/2 + 3/2)' 

where, as in Eq. (f3"5)) . we assumed that q e oc £)"" d / 2 ; 
and 6pp{ (k, 77) = g a b{ii) </>b(k) is the usual linear the- 
ory evolved density field in the PPF approximation. 
We can then write the density power spectrum to 
leading order in PPF corrections as 



P ss {k) = Pppf(k)- 



(rw/2-l)(rw/2 + 3/2) 



, (45) 
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where P pp f(k) is the linear density power spectrum 
in the PPF approximation, and Ps qa (k) is the cross 
power spectrum given by 

(5(k) ge (q)) = 5 D (k+q) P, w (k), (46) 

which we measure from the numerical simulations. 
From Eq. (|42p we can also obtain the velocity diver- 
gence in linear theory, 



Qe(k, rj) (n vd /2) 



(n vd /2- l)(n„ d /2 + 3/2)' 
(47) 

and the corresponding power spectrum, 



P ee (k) = P PP f(k)- 



f 9 g„ (k) 



(n vd /2-l)(n vd /2 + 3/2)' 



(48) 



Note that the correction in the case of the velocity 
divergence power spectrum is a factor of (n vd /2 m 3) 
larger than in the case of the density; one can under- 
stand this from the fact that the divergence responds 
to the rate of change of the density fluctuations, and 
the correction to S pp { grows as D™ 1 "^ 2 . It is also 
worth noticing that the corrections are negative, i.e. 
velocity dispersion tends to reduce the growth of 
structure; this is also expected for a stress tensor 
paij w —pSij with a positive pressure (due to ther- 
mal motions) that is positively correlated with den- 
sity fluctuations. Figure [5] shows the absolute value 
of these corrections relative to the PPF approxima- 
tion for both density and velocity divergence at red- 
shifts z = 0, 0.5, 1. We see that at z = the correc- 
tion to the divergence power spectrum reaches 1% at 
about k ~ 0.1 ZiMpc , while for the density power 
spectrum this happens at about k ~ 0.2 /i Mpc -1 . 
By z = 1 these scales shift by about a factor of 
two. At higher redshifts they rapidly decline as the 
growth factor changes rapidly before the onset of 
cosmic acceleration. 

These effects can be understood qualitatively, and 
to some extent quantitatively, by considering the 
typical size of the corrections to the velocities pre- 
dicted by the single-stream, PPF approximation, 
smoothed over scale R. These corrections are, in 
average, of order 

1/4 



(49) 



a r v ms (R)=nf d 6 kP a {k)W^{kR) 




FIG. 10: Top panel: root mean square position fluctua- 
tions, Eq. (|50p . induced by velocity dispersion smoothed 
at scale R divided by R. Bottom panel: rms veloc- 
ity dispersion, Eq. (|49[1 . in solid lines compared to rms 
bulk motions (dashed lines) smoothed on scale R. Note 
that velocity dispersion is smoothed on scales of order 
1 ft" 1 Mpc, thus the solid line is an underestimate at 
small scales. All the quantities in this figure are eval- 
uated at z = 0. 



where Wth(^-R) is the Fourier transform of a top- 
hat filter of radius R, P a {k) is the power spectrum 
of the trace of the velocity dispersion tensor (which 
is the dominant component), and the factor Ti.f re- 
stores the correct units to <Ty (see Eq. IMl) . Equiva- 
lently, these velocity corrections can be interpreted 
as comoving position fluctuations, given by 



d 3 fcP CT (fc) Wf H (feR) 



1/4 



(50) 



These two quantities are shown in Fig. [TO] In the 
top panel, the ratio of the displacement corrections 
from Eq. ([50]) to the scale R is plotted as a function 
of scale. An order of magnitude estimate of the effect 
on the density power spectrum can be obtained from 
the following argument. The dispersion in comoving 
positions given by A rms (R) smooths out density per- 
turbations. That suppression is approximately given 
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by 

P sm ooth(fc) ~ P(k)e'^ kArms ^/ k ^ 2 . (51) 

At large scales, e.g. k ~ 0.1 /iMpc -1 , this gives 
a suppression consistent with the previously calcu- 
lated density power spectrum corrections seen in 
Fig. M 

The bottom panel shows a" as as a function of 
R, Eq. (|49|) . We can see that the velocity disper- 
sion on scales of ~ 100 h^ 1 Mpc is of order 15 km/s. 
Comparing this dispersion with the single-stream 
bulk velocities on the same scale (dashed line), we 
conclude that the velocity dispersion corrections on 
those scales are small but, nevertheless, larger in rel- 
ative terms than for the density power spectrum, in 
agreement with the detailed calculation presented in 
Fig. US 

In [551 ] it was argued that percent level corrections 
from orbit crossing to the density power spectrum 
are expected at k ~ 0.1 ft, Mpc -1 based on a model 



of "sticky dark matter" . The effect discussed in that 
work is not an estimate of deviations from the PPF 
approximation, see Appendix A in [ii^ for more de- 
tails. Here we note that the estimate in [55[ for 
the density power spectrum is two times larger than 
found here, and opposite in sign. 



C. PT + Vorticity 

As discussed above, the effects of the vector modes 
of the stress tensor (q w ) on the density and di- 
vergence power spectrum only appear beyond lin- 
ear theory, since in linear theory scalar and vector 
modes are decoupled. Here we estimate these correc- 
tions by calculating the one-loop density power spec- 
trum including the effects from the vorticity, which 
is sourced by q„,. The strategy is as follows. We 
rewrite the equations of motion for density pertur- 
bations including now the nonlinear terms as follows, 



d n 9(k) 



d v 6(k) - 0(k) = / d 3 hd 3 k 2 S D (k - k 12 ) 



0(k) 3<J(k) 



d 3 k 1 d 3 k 2 S D (k-k 12 ) 



k • k 2 ki x k 2 

-T2-M2 72— ■ *1W2 



fc 2 (k! • k 2 )M 2 (ki • k 2 )(ki x k 2 ) • # lW2 



91,21,2 



h 2 h 2 



(52) 
(53) 



r 



where ki 2 = ki + k 2 . To avoid cumbersome expres- 
sions, we have not written the time dependence of 
the fields explicitly. Also, on the right hand side, the 
subscripts "1" and "2" mean the fields evaluated at 
ki and k 2 respectively. 

A couple of points are worth noting. We have de- 
composed the velocity field into a divergence and a 
vorticity, and the latter will be taken as a known 
forcing term in the equations, since w can be solved 
in linear theory as an uncoupled field from the mea- 



sured (\ w (Eq. I37[) or directly measured from the 
Delaunay method. In addition, note that we have 
neglected the qg source in the equation of motion 
for 9, since this effect was included already in the 
last section; here we are only interested in correc- 
tions due to q w alone, which enter through the w 
forcing terms. 

Following the compact notation introduced in the 
previous section, we can rewrite Eqs. (|52H53j) as 



d v ip a (k,ri) +Q ab il} b (k,r)) = J d 3 kid 3 k 2 S D (k - ki 2 )[7 ahc (ki, k 2 )?/>&(ki, 7y)?A c (k 2 , 77) + A ab (ki, k 2 , 77)^6^1, r))], 

(54) 
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where 7 characterizes the non-linear mode coupling amplitudes and A is the w-dependent forcing term. 
They can be written as 



7112 



ki2 ■ k 2 



7222 



fci 2 (ki • k 2 



ki x k 2 

^ab = 73 w(k 2 ,?7) 



1 
n k i k 2 

u -sr 



(55) 



Equation (|54[) can be formally solved by introducing again the linear propagator g a b, yielding: 
^a(k, 77) = g a b(r))(f>bQs.) + ( dsg ab (r)-s) f d 3 kid 3 k 2 S D (k - k X2 ) b^k]., k 2 )V' c (k 1 , s)V>d(k 2 , s) 



+ J 4b c (ki,k 2) s)V'c(ki,s)] 



(56) 



This is an implicit form of the solution - it is written in terms of itself. However, it is suitable for a 
pcrturbative method. We write the field tjj as a perturbative series: 



(57) 



and combining Eqs. (|5S|) and ([57]) . we get a solution for the n-th order fields in terms of the lower order 
fields: 



^i n) (M) 



dsga^-s) d 3 k x d 3 k 2 (SD(k-k 12 )k C(1 (ki,k 2 ) Yl ^ r) (s> k i)^ S) (s.k 2 ) 



r+s— n 



+A fec (k 1 ,k 2)S )^"- 1 )(k 1)S ) 



(58) 



Thus, by knowing the linear solutions, we can calculate the solutions to any order. The linear solutions for 
S and 9 are just the PPF linear solutions. In order to solve for the higher order fields 5 and 6, we need the 
vorticity field and its time dependence, which we have measured from dark matter N-body simulations. 

As discussed above, the leading order correction due to vorticity effects appears in the one-loop contribution 
to the power spectrum. The density power spectrum, to that order, can be written as 



P S (k)5 D (k + q) = (^(k^^q)) + [<^ 2 >(k)^ 2 >(q)) + 2<<5« (k)<^ (q)}* 



(59) 



The first term is the usual tree-level (linear theory) power spectrum, while the terms in square brackets 
correspond to the one-loop correction, which are usually written as P 22 + Pi 3 due to their dependence on PT 
order. Let us start by focusing on the first one-loop term, which leads to the P 22 contribution. The second 
order density can be written as 



(60) 



<5< 2 >(k) = / d 3 q 1 d 3 q 2 5 D (k-< ll2 ) [P(qi,q 2 )<5( 1 Hqi)5 (1) (q 2 ) + G(qi,q 2 )-w(q 1 )^ 1 )(q 2 ) 
= J d 3 qid 3 q 2 5 D (k - q i2 ) [Fi 2 <5i<5 2 + G\2 ■ wi<5 2 ]. 



Note that in the second line we have just simplified the notation. The P 22 contribution to the power spectrum 
has four terms: 



(61) 



P 22 (fc)M k + q) = J d 3 q 1 ...d 3 q i S D (k-q 12 )S D (ci-ci 3i )^F 12 F 3 4 : {5id 2 d 3 S i ) 

+P12G34 • (M 2 W3<5 4 ) + G 12 P 34 • (wi&W + G^G^Kc^f^} 
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The first term is the usual PPF P 22 one-loop contribution. The second and third terms contain factors of 
the form (<5w), which vanish due to symmetry. Then at large scales where connected contributions can be 
neglected we have 

Kfctflffc) = (w?v4 )(S 2 S 4 ) = MQi + q 3 ) M<fc + q 4 ) P^(qi)P 5 (q 2 ) (62) 
from which we get the vorticity contribution to the P 22 power spectrum: 

AP 22 (k) = J AG^-q,q-k)G Q (q,k-q)P^(q)P 5 (k-q) = J #q \ G{r ** ~ ^ P w {q)P s (\k - q|), 

(63) 

where we have used the actual vector structure of G(k, q) and the fact that 



pa/3, \ = p w{o) 

WW \H/ rt 



o a p 5— 

T 



(64) 



Similarly, we can compute the Pi 3 contribution to the density power spectrum. The third order density 
field can be written schematically as 

<5 (3) (k) = J d i q l d z q 2 d z q^5 D ^ - ci 12Z )[H{c{ l ,ci 2 ,ci Z )5 1 825? > + R{q.i,^ 

(65) 

Then, the vorticity contribution to the power spectrum reads 

AP 13 (fc) = P 5 (k) J Ax^(qi,-qi,-q)^(qi) = Ps{k) J d 3 q S aa (q, -k)P w (q), (66) 

where in the last equality we have used Eq. . If one assumes that the time dependence of the vorticity 

is given by w a D™"^ 2 , as found in section [Til it is possible to write explicit expressions for the AP 22 and 
AP13 power spectra: 

AP m /inMpfl. h 2fc 2 (l-z 2 ) [(3 + n w )k 2 + (1 + n w )q 2 - 4(1 + n w /2)k ■ q] 2 
AP 22 (fc) = J d qP w ( q )Ps(\k - q|) - 2 X n2u(5 + nw)2 | k - q | 2 • ( 6? ) 



AP 13W = -Ps(k) j <P qPw{q) g2|k + J^J ){5 + 2nw) {* [0 + n w )(3 + 2n,)*» 

+ (3 + 9n w + 2n 2 w )q 2 } + 2(1 + n w ) k ■ q [q 2 + (9 + 2n,„)fc 2 ] + 4(2 + 3n TO /2)(k • q) 2 }, (68) 
where x is the cosine of the angle between k and q. 



The end result of these calculations is that the 
leading large-scale contribution of vector modes of 
the stress tensor to the density power spectrum is 
fully specified in terms of the autocorrelation or 
power spectrum of the vorticity, which we have mea- 
sured from the simulations. Figure [11] shows the 



results of these calculations. We see that the to- 
tal correction is negative, as expected physically, 
and very small at large scales. For example, at 
fc ~ 0.1 h Mpc -1 where the scalar modes contributed 
percent level corrections, the modifications of the 
PPF approximation from vector modes is about 
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FIG. 11: Corrections to the density power spectrum at 
z — due to stress tensor vector modes (vorticity ef- 
fects), see Eqs. (J67J and ((68j. Note that the AP i3 contri- 
bution (long dashed lines) is negative and larger in mag- 
nitude than the AP22 contribution (dashed lines). The 
total correction (solid lines) is negative and reaches 1% of 
the linear spectrum (top dotted lines) at k ~ 1 AMpc , 
where further nonlinear effects not included here should 
become important. 



10~ 4 of the linear spectrum, and thus totally neg- 
ligible. The reason for this is that by symmetry 
the vorticity does not couple to the scalar modes, 
it is only through vorticity squared that the effect is 
present. We expect similar results for the velocity 
divergence power spectrum within a factor of a few, 
still completely negligible at large scales. 



V. CONCLUSIONS 

We studied the impact of orbit crossing in the 
large-scale power spectra of density and velocity di- 
vergence fields, which are usually described in the 
pressureless perfect fluid (PPF) approximation. We 
presented a method to extend perturbation theory 
(PT) beyond the PPF approximation, based on mea- 
suring the stress tensor induced by orbit crossing in 



numerical simulations. The stress tensor, when de- 
composed into scalar and vector modes leads to cor- 
rections associated with velocity dispersion and the 
effects of vorticity. We found the effects due to the 
scalar modes to be small, but not negligible at large 
scales (k ~ 0.1 /iMpc -1 ), particularly for the ve- 
locity divergence power spectrum (see Fig. The 
impact of vorticity on large scales is much smaller, 
see Fig. QT] These two effects appear at different or- 
ders in PT and have been included separately as we 
are interested in large scales where the induced cor- 
rections are small. Both lead to suppressions of the 
power spectra predicted by the PPF approximation, 
as expected physically since velocity dispersion and 
vorticity should inhibit collapse. In this regard we 
emphasize that neglecting orbit crossing has oppo- 
site effects on Eulerian compared to Lagrangian PT. 
For Lagrangian PT, neglecting orbit crossing leads 
to (much more severe) underestimates of the density 
power spectrum (see e.g. [1] for a recent example), 
since neglecting self-gravity in caustics leads to arti- 
ficial thickening of such structures when trajectories 
cross without interacting. 

A novel aspect of our calculation is the estimation 
of the stress tensor and the vorticity and divergence 
power spectra from numerical simulations. To esti- 
mate velocity fields, we applied the Delaunay tessel- 
lation method, which we have shown to be a more 
reliable estimator than traditional mass weighting 
schemes. While estimates of the velocity divergence 
are robust, we found that measurements of the vor- 
ticity power spectrum are significantly more diffi- 
cult, due to aliasing during the measurement pro- 
cess and most importantly lack of resolution in the 
simulations. For the latter we have found that low 
resolution simulations can overestimate the vortic- 
ity power spectrum by an order of magnitude. This 
maybe be due to insufficient spatial resolution in 
multistrcaming regions, with the overestimate per- 
haps related to aliasing effects during the PM part 
of the force calculation, which may generate a vec- 
tor mode. In any event, for high enough resolution 
we find that the vorticity power spectrum converges 
to a stable answer. On the other hand, care must 
be taken that these spurious effects are not present 
when using numerical simulations to study nonlin- 
ear velocities, since artificial vorticity can amplify 
the velocity power spectrum at small scales. 

A nontrivial check of our numerical calculation 
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of the stress tensor, which we have done using an 
adaptive method independent of the Delaunay tes- 
selation, is that its vector modes source the growth 
of vorticity. Therefore, using linear PT from this 
vector source one should recover at large scales the 
vorticity power spectrum measured by the Delaunay 
method, as we do (see Fig. [HJ) . This does not test the 
scalar mode of the stress tensor though, which ends 
up inducing the largest correction to the PPF ap- 
proximation. In this respect, it would be interesting 
to test how robust the scalar part of the stress tensor 
is to details of the numerical simulations, as spuri- 
ous effects due to discreteness may amplify veloc- 
ity dispersion in simulations [56|, [HtJ (see also (58j). 
As far as we know, our work is the first to make a 
quantitative connection between the growth of ve- 
locity dispersion and that of the density power spec- 
trum, which will be useful to probe more in order 
to make sure that simulations can correctly repro- 
duce the matter power spectrum to percent level, 
as required for the next generation of weak lensing 
surveys designed to probe cosmic acceleration [59(. 

The deviations we found from the PPF approxi- 
mation at large scales are small but not negligible, 
in particular for the velocity divergence power spec- 
trum, for which corrections are a factor of about 
three larger than for the density power spectrum. 
Our estimate, being based on numerical simulations, 
corresponds to fixed cosmological parameters (e.g. 
(Ts = 0.9, Q m = 0.27 and n s — 1). Given the strong 
dependence on the growth factor of the correction 
(oc D 2 : 25 relative to PPF) we expect it to be smaller 
for lower normalization amplitudes, as well for cos- 
mological parameters that correspond to less power 
at small scales (i.e. il m < 0.27 and n s < 1). 

In section UTTl we sketched what must be done to 
include these effects from first principles into ana- 
lytic calculations such as RPT that usually start 
from the PPF approximation, instead of using the 
hybrid approach we develop here partially based on 
numerical simulations. Including velocity dispersion 
in a self-consistent manner should cure divergences 
that appear in PT for scale-free models with initial 
power spectra P(k) ~ k n for n > —1, as well as 
regulate the divergences that appear in the resum- 
mation of the Lagrangian space propagator [6(| ■ In 
addition, we expect velocity dispersion and vorticity 
to be crucial to describe the virial turnover in the 
density power spectrum. Another interesting appli- 



cation of the ideas presented in section IIIII is to use 
the cumulant hierarchy to describe nonlinear effects 
in a massive neutrino component, to improve on re- 
cent calculations [6l|, |62| that assume linearity. We 
hope to report on some of this in the near future. 
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APPENDIX A: ESTIMATING THE 
VELOCITY FIELD 

1. The Delaunay Tessellation 

One of the main obstacles in measuring the ve- 
locity field from cosmological simulations is the fact 
that it is only sampled on a discrete set of points. 
One encounters the same difficulty when recon- 
structing the peculiar velocity field from observa- 
tions, where velocity is only sampled at the locations 
of galaxies. One can identify two problems associ- 
ated with this fact. On one hand, since the velocity 
is only known at points where the mass is located, al- 
most all procedures to reconstruct the velocity field 
from a discrete sample give mass-weighted quan- 
tities, while most theoretical predictions concern 
volume-weighted quantities. On the other hand, low 
density regions are very sparsely sampled, and there- 
fore subject to large Poisson errors. Laying out a 
structured grid, as it is often done, leaves grid points 
empty in such regions with the consequence that 
the velocity field is set to zero (thus missing out- 
flows in voids), while in practice the velocity field 
is undetermined due to the poor mass resolution. 
These issues have been recognized for a long time in 
the theoretical large-scale structure literature, see 
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e.g. isai2ai2si- 

The work in [37j introduced two new velocity es- 
timation methods that attempt to overcome these 
problems. These methods are based on the Voronoi 
and Delaunay tessellations of the discrete set of 
points where the velocity field is sampled. They 
showed that the Delaunay Tessellation method has 
fewer computation requirements than the Voronoi 
Tessellation method, while giving equally or more 
reliable results. Thus, in our work we will consider 
only the Delaunay method, following [37] to a large 
extent. Our implementation of the method is based 
on the public code Qhull [40J to construct the De- 
launay tessellation. Many other applications of De- 
launay tessellations have recently appeared in the 
large-scale structure literature (see e.g. [HI, EH), see 
also [43j for an in-depth review and other applica- 
tions. 

The formal definition of the Delaunay tessellation 
T>{V) of a set of points V (in three dimensions) is 
the set of tetrahedrons defined by four points whose 
circumscribing sphere is empty in the sense that no 
point of the generating set V should be inside the 
circumsphere In Fig. [T2] we show an example 

of the tessellation of a set of random points in two 
dimensions. It can be shown that the Delaunay tes- 
sellation is unique. Moreover, the Delaunay tetra- 
hedrons are objects of minimal size and elongation. 
These characteristics make the Delaunay method op- 
timal for a three dimensional interpolation. 



2. Reconstructing the Velocity Field 

Once the Delaunay tessellation from the set of 
sample points is obtained, it is possible to estimate 
the velocity at each point p in space by linearly inter- 
polating the velocities at the vertices of the tetrahe- 
dron that contains the point p. This procedure leads 
to a continuous velocity field with constant gradient 
in each tetrahedron. 

Mathematically, we can express the Delaunay 
method to find the velocity u at a point p with co- 
ordinates x as follows. Let Xj, with i = 0,1,2,3, 
be the coordinates of the vertices of the tetrahedron 
containing p. Since the Delaunay tetrahedrons are 
non-degenerate (i.e. they do not collapse into 2D 
objects), we can express x as a linear combination 



of x,-: 



Ax = «jAxj 



(Al) 



,.=i 



where Ax = x — xo and Ax^ = x, — xo. The linear 
interpolation of the velocity at point p is simply 



Au = y^ajAui, 



(A2) 



;.=i 



where Au = u — Uo, Au, = — uo, and on satisfy 
Eq. (|A1[) . Thus, the problem reduces to solve for the 
cti, which can be readily be written as 



(A3) 



where Ax, Ay and Az are the components of x, and 
A consists of the components of the Ax^: 







f Ax \ 


a 2 




Ay 


\a 3 J 




{AzJ 



( Axi Ax 2 Ax 3 
Ayi Ay 2 Ay 3 
\ Azi Az 2 Az 3 



(A4) 



These equations allow us to compute an estima- 
tion of the velocity field for any point in the simula- 
tion volume. In particular, we are interested in de- 
termining the volume-averaged field at a given set of 
points, often a grid. Typically, one wants to compute 
the average u^(r^) of the velocity field in spheres of 
radius R centered at the points r^, usually on a grid. 
In order to obtain that, one can carry out the fol- 
lowing algorithm [37| : 

1. Construct the Delaunay tessellation from the 
locations of the simulation particles. 

2. For each point r^: 

(a) Find the intersection of the Delaunay 
tetrahedrons with a sphere of radius R 
centered at r^. 

(b) For each intersecting tetrahedrons j, de- 
termine intersection volume Vj and mean 
velocity Uj in that volume. 

(c) Compute V VjUj /(AttR 3 /3). 
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FIG. 12: Delaunay tessellation of a set of random points in two dimensions. The left panel shows the original set of 
points. In the right panel, we show the regions corresponding to the tessellation. Note that in two dimensions the 
Delaunay tessellation consists of triangles instead of tetrahedrons. 



However, both constructing the tessellation for the 
large number of particles (~ 10 9 ) typical of state-of- 
the-art simulations and finding the tetrahedrons in- 
tersecting a given sphere are very time-consuming. 
For the sake of efficiency, we modified the previous 
procedure. Instead of calculating the velocity aver- 
age on a sphere centered at a grid point rj , we com- 
pute the average in the volume given by all Delaunay 
tetrahedrons totally contained in the grid cell corre- 
sponding to Yi. Thus, our approximate algorithm 
reads as follows: 

1. For each grid point r*, 

(a) Construct the Delaunay tessellation of 
the points contained in the corresponding 
grid cell. 

(b) Compute the volume Vj and mean veloc- 
ity Uj for every Delaunay tetrahedron. 

(c) Compute £j VjUj/Y^j Vj. 

It is possible to write explicit expressions for the vol- 
ume Vj and mean velocities Uj in a Delaunay tetra- 
hedron. It follows from elementary geometry that: 

V$ = |det(A)|, u^i^Tuf, (A5) 

4=0 



where A is defined in Eq. (|A4[) , and uv are the four 
velocities at the vertices of the tetrahedron j. 

Note that in the new algorithm, we only construct 
the tessellation of the points inside the grid cell, a 
much smaller number of particles than the total sim- 
ulation. Moreover, it is no longer necessary to find 
the tetrahedrons or their intersection with a sphere. 
Nevertheless, with this procedure, the average vol- 
ume will in general vary from grid point to grid 
point. We reduce this undesired effect by only using 
relatively coarse grids, where we expect a more uni- 
form distribution of particles. We are thus obtaining 
a smoothed field ur(x) given by 

Ufl(x) = J d 3 y W R {x - y) u(y), (A6) 

where Wr(x) is a spherical top- hat filter with R w 
(?>V CC \\/ kit) 1 / 3, . To deconvolve Fourier space quan- 
tities, we thus divide by the Fourier transform of 
Wr. Note that this correction is only correct on av- 
erage, tests reveal that it is accurate to about 1% at 
k = 0.2 h Mpc -1 , more than enough for our purposes 
in this paper, but not enough for precision tests of 
velocity divergence power spectra in the weakly non- 
linear regime. 
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FIG. 13: Left panel: Delaunay estimated divergence and vorticity power spectra. The solid lines correspond to the 
power spectra used to generate the velocity field. The dashed lines are the spectra measured on a 120 3 grid, and the 
dotted lines are measured on a 76 3 grid. All spectra are corrected by deconvolving the smoothing kernel (Eq. IA6|l . 
Right panel: Same as left panel, but using the CIC method. Note that since these measured spectra are obtained as 
the ratio of two interpolated quantities, they cannot be easily corrected for the window of the interpolation scheme. 



3. Testing the Delaunay Method 

One of the difficulties of testing accuracy of the 
Delaunay method is that it is expected to be more 
accurate (in measuring volume- weighted quantities) 
than the traditional estimations from mass-weighted 
schemes. Thus, we lack a more trustworthy method 
to use as reference. Out strategy to overcome this 
difficulty is setting up a "controlled numerical ex- 
periment" : we generate a random Gaussian velocity 
field with given divergence and vorticity power spec- 
tra and then use the Delaunay method to recover the 
velocity statistics. In this way, we can compare the 
results of the method with the exact input power 
spectra used to generate the velocity field. 

For the sake of comparison, we also measure the 
velocity power spectra with the well known Cloud- 
in-Cell mass- weighted method (CIC). It consists of 
interpolating the particles mass and velocity on a 
grid using the CIC kernel Wcic(x) = Hi Wcic(xi) 



defined by 

W cl c( Xi ) = \ 1 - ]Xil ft*';? (A7) 
{ for \xi\ > 1, 

where x is measured in units of grid separation. Note 
that interpolating the particle velocities by using 
this method gives the momentum field instead of 
the velocity field. Thus, one needs to compute the 
ratio between this quantity and the density field to 
obtain the velocity. As we mentioned above, in un- 
derdense regions if the grid is made too fine there 
will be grid points for which no particle is assigned, 
which means there is no information on the velocity 
field, but typically one would set to zero (incorrectly) 
the velocity. In addition, dividing the interpolated 
momentum by the interpolated density means that 
it is difficult to correct for the interpolation kernel 
after the velocity field is Fourier transformed, unlike 
the case of the density field (see [45[ for a recent dis- 
cussion of interpolation corrections and comparison 
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of CIC with other mass assignment schemes for the 
density field). 

We generate a Gaussian velocity field on a grid 
of 400 3 cells with a divergence and vorticity power 
spectra based roughly on expectations from previous 
measurements in the literature [46[ and then inter- 
polate this velocity on the positions of the 640 3 dark 
matter particles obtained from running an N-body 
simulation with Gadget2 (4?| (see Table fl] below for 
more details on the simulations). Then we measure 
the divergence and vorticity power spectra using the 
Delaunay method and the CIC method. 

The results are shown in Fig. [T21 We applied the 
Delaunay method, as described in the previous sub- 
section, on a coarse grid of 76 3 cells and a finer grid 
of 120 3 cells, and measured the power spectra us- 
ing fast Fourier transforms. On scales close to the 
Nyquist frequency, the power spectra were corrected 
by deconvolving the kernel defined in Eq. (|A6|1 . The 
recovered divergence and vorticity agree very well 
with the input power-spectra (left panel). In the 
right panel, we show the results of the CIC method. 
We see that in order to obtain results comparable 
to the Delaunay method, one needs to use a much 
finer interpolation grid. Even in that case, there are 
considerable differences on large scales. On small 
scales, the power spectrum is underestimated. This 
is due to the CIC interpolation kernel. However, 
it cannot be corrected as in the Delaunay case be- 
cause the velocity field was obtained as the ratio of 
two CIC-interpolated quantities. In principle, one 
could deconvolve the density and the momentum in- 
terpolated fields before taking the ratio, but this pro- 
cedure does not give good results because it intro- 
duces noise in the deconvolved fields 26] . In partic- 
ular, we observed that the deconvolved density field 
has a non-negligible number of negative density grid 
points. Similar results on the better noise properties 
of the Delaunay method were obtained for the PDF 
of the velocity divergence in (37| . 

4. Sampling Effects: Aliasing 

One of our goals in this paper is to estimate re- 
liably the vorticity field from cosmological simula- 
tions. Since at large scales it is expected to be very 
small compared to the divergence, it is important to 
analyze the sampling effects. We will show that such 



effects make the sampled vorticity field a mixture of 
the vorticity and divergence of the original field. 

Let us assume we know the velocity field u(x) in- 
side a box of volume L 3 , and will study the effects 
of sampling that field on a grid of N 3 cells. Let us 
decompose the original (exact) field u(x) in Fourier 
series: 

u(x) = ^u(k)e lkx , (A8) 

k 

where the vector k has components k a = 2irn a /L 
with n a 6 Z, i.e. arbitrarily large frequencies appear 
in the Fourier sum. We want to compare these exact 
Fourier modes u(k) to the discrete Fourier modes 
u(q) on the grid (q a — 2irm a /L, m a £ Z, < m a < 
N). We can write the discrete Fourier transform of 
the velocity field as 

x 

= ^EE a ( k ) e * (k ~ q)x < ( A9 ) 

x k 

where q and x are on the grid. This can be further 
simplified into 

The product vanishes unless r = k— q = with 
m an integer vector. Finally, we obtain 

u (q) = z^ a ( c t + r )' r = - ^r m - ( An ) 

r 

This equation states explicitly that velocity Fourier 
modes beyond the Nyquist wavenumber of the grid 
(fcjvj, = irN/L) affect the grid Fourier modes. This 
is known as aliasing. To see how this effect appears 
in the vorticity power spectrum, let us assume that 
the velocity field is purely potential, that is, u(k) = 
ik9(k)/k 2 . Evidently, w(k) = ik x u(k) = 0, but 

i(q)Eiqxu(q) = ^^(q + r) (A12) 

does not vanish. Moreover, the sampled vorticity 
power spectrum can be written as 

^W=£jf^jl P o(h + r\), (A13) 
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FIG. 14: Measured divergence and vorticity power spec- 
tra from a vorticity-free velocity field sampled on three 
different grids of 20 3 , 50 3 and 100 3 grid points (top 
to bottom). This shows that the measured vorticity is 
purely due to aliasing. The dotted lines correspond to 
the predictions of Eq. (| A13f) . 



where r = (2nN/L)m denotes multiples of the 
Nyquist frequency. Note that Pg is the power spec- 
trum of the original divergence field. 

Equation (IA13[) tells us that the velocity diver- 
gence of wavenumbers larger than the Nyquist of 
the grid induces a spurious vorticity in the sam- 
pled field. In the low-q limit (q <C 2ttN/L), it is 
easy to see that P w (q) oc q 2 . Figure [14] shows the 
predictions of the large scale limit of this formula. 
We generated a zero-vorticity velocity field, which 
was sampled without smoothing it on different grids. 
Then, we compared the measured vorticity power 
spectrum to the predictions of Eq. (|A13[) . Note that 
the amplitude of the correction was not fitted: the 
formula estimates correctly both the low-/c limit de- 
pendence of the vorticity power spectrum and the 
amplitude of the spurious vorticity. 

It is important to remark that we cannot directly 
extrapolate these results to the vorticity estimates 
of Fig. [T2J In that case, both the vorticity and 
divergence fields are smoothed on a scale given by 
the grid separation. This procedure greatly reduces 
the power spectrum on wavenumbers larger than the 
Nyquist frequency of the grid, making the aliasing 
effect less important. However, this analysis is useful 
to set an upper-bound estimate of aliasing effects. 
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